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Using 7 -ray data from the Fermi Large Area Telescope, various groups have identified a clear 
excess emission in the Inner Galaxy, at energies around a few GeV. This excess resembles remarkably 
well a signal from dark-matter annihilation. One of the most compelling astrophysical interpretations 
is that the excess is caused by the combined effect of a previously undetected population of dim 
7 -ray sources. Because of their spectral similarity, the best candidates are millisecond pulsars. 
Here, we search for this hypothetical source population, using a novel approach based on wavelet 
decomposition of the 7 -ray sky and the statistics of Gaussian random fields. Using almost seven years 
of Fermi-UKT data, we detect a clustering of photons as predicted for the hypothetical population 
of millisecond pulsar, with a statistical significance of 10.Oa. For plausible values of the luminosity 
function, this population explains 100% of the observed excess emission. We argue that other 
extragalactic or Galactic sources, a mismodeling of Galactic diffuse emission, or the thick-disk 
population of pulsars are unlikely to account for this observation. 


Introduction. Since its launch in 2008, the Fermi 
Large Area Telescope (LAT) has revolutionized our un¬ 
derstanding of the 7 -ray sky. Among the major successes 
are the detection of more than 3000 y-ray sources [ij, 
the discovery of the Fermi bubbles Q , some of the most 
stringent limits on dark-matter annihilation Q and, most 
recently, the detection of cross-correlations between the 
extragalactic 7 -ray background and various galaxy cata¬ 
logs y. 


One of the most interesting 7 -ray signatures identified 
in the Fermi-LPGT data by various groups HE! , is an 
excess emission in the Inner Galaxy at energies around a 
few GeV. This excess attracted great attention because 
it has properties typical for a dark-matter annihilation 
signal. This Galactic center excess (GGE) is detected 
both within the inner 10 arcmin of the Galactic Center 
(GC ) H , ^ 13 and up to Galactic latitudes of more than 
10 ° HillinBlil- It features a remarkably uniform 
spectrum and approximately spherical symmetry 13 . .[i3- 
Proposed diffuse emission mechanisms, like leptonic or 
hadronic outbursts |19l42ll or cosmic-ray injection in the 
central molecular zone [ 22 | . potentially explain part of 
the excess emission. However, it is challenging to explain 
all of the above aspects of the GGE with these mecha¬ 
nisms alone. 


Probably the most plausible astrophysical interpreta¬ 
tion for the GGE is the combined emission from a large 
number of unresolved millisecond puls ars (MSPs) in the 
Galactic bulge region [13 E, E,Hi. MSPs feature a 
spectrum compatible with the GGE emission [l^, and 
a large unresolved component can naturally explain the 
uniformity of the GGE spectrum in different regions of 
the sky. Recently, it was shown that the spatial distribu¬ 
tion of MSPs that were spilled out of disrupted globular 
clusters can explain the morphology of the GGE [25| . 
Such MSPs from disrupted globular clusters have also 
been suggested as the source behind the GeV through 
TeV emission in the inner few parsec of the GC (E] . Fur¬ 


ther possible support for the MSP hypothesis might come 
from Chandra observations of low-mass x-ray binaries 
(which are progenitor systems of MSPs) in M31, which 
show a centrally peaked profile in the inner 2 kpc 27, 
as well as the recent observation of extended hard X-ray 
emission from the Galactic Center by NuSTAR (^ . 

It was claimed that an interpretation of 100% of 
the GGE emission in terms of MSPs would be already 
ruled out: a sizeable fraction of the required lO^-lO"* 
MSPs should have been already detected by the Fermi- 
LAT [13 El, but no (isolated) MSP has been identified 
so far in the bulge region. This conclusion depends cru¬ 
cially, however, on the adopted 7 -ray luminosity of the 
brightest MSPs in the bulge population, on the effective 
source sensitivity of Fermi-GPUP, and on the treatment of 
unassociated sources in the Inner Galaxy miai. A real¬ 
istic sensitivity study for MSPs in the context of the GeV 
excess, taking into account all these effects, was lacking 
in the literature up to now (but see Ref. [SSj)- 

In this Letter, we close this gap and present a novel 
technique for the analysis of dim 7 -ray sources and ap¬ 
ply it to Fermi-PiPGP observations of the Inner Galaxy. 
Our method is based on the statistics of maxima in the 
wavelet-transformed 7 -ray sky (in context of Fermi-PPUP 
data, wavelet transforms were used previously for the 


a previor 

identification of point source seeds [l|, [la])- We search 
for contributions from a large number of dim MSP-like 
sources, assuming that they are spatially distributed as 
suggested by GGE observations. Our method has sev¬ 
eral advantages with respect to previou^ proposed tech¬ 
niques based on one-point fluctuations |34| , most notably 
the independence from Galactic diffuse emission models 
and the ability for candidate source localization. 


Modeling. We simulate a population of MSP-like 
sources, which we hereafter refer to simply as the cen¬ 
tral source population (CSP), distributed around the GC 
at 8.5 kpc distance from the Sun. The CSP is taken to 
have a spatial distribution that follows a radial power law 
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with an index of F = —2.5 and a hard cutoff at radius 
r = 3 kpc [13, [I3|- As a reference 7 -ray energy spec¬ 
trum, we adopt the stacked MSP spectrum from Ref. [35| . 
^ oc e-'®/3-78GeV^-i.57^ luminosity func¬ 


tion is modeled with a power law, ^ oc L~°‘, with index 
a = —1.5 32l. l35l - l37l| . and with lower and upper hard cut¬ 
offs at Fmin = 10 ^® ergs“^ and Fmax = 10 ®^- 10 ®®ergs“^, 
respectively. Luminosities are integrated over 0.1-100 
GeV. Our results depend little on Lmin- Given that 
only about 70 MSPs have been detected in 7 rays up 
to now [ 33 , Lmax is not well constrained. The 7 -ray lu¬ 
minosity of the brightest observed MSP is somewhere in 
the range (0.5-2) x 10®® ergs“^ depending on 

the adopted source distance 2^ l32| . Diffuse emission is 
modeled with the standard model for point source ana¬ 
lysis gll_iem_v06. f its and the corresponding isotropic 
background. 


Data. For our analysis, we use almost seven years of 
ultraclean Fermi-LAT P8R2 data taken between August 
4 2008 and June 3 2015 (we find similar results for source 
class data). We select both front- and back converted 
events in the energy range 1-4 GeV, which covers the 
peak of the GCE spectrum. The region of interest (ROl) 
covers the Inner Galaxy and spans Galactic longitudes 
\(.\ < 12° and latitudes 2° < | 6 | < 12°. The data are 
binned in Cartesian coordinates with a pixel size of 0.1°. 


Wavelet peaks. The wavelet transform of the 7 -ray 
data is defined as the convolution of the photon count 
map, C(D), with the wavelet kernel, W(D), 


Dw[c]in) = j dn'Win - n')c{n'), (i) 

where D denotes Galactic coordinates [note that 
f dnW(n) = 0]. The central observable for the current 
analysis is the signal-to-noise ratio (SNR) of the wavelet 
transform, which we define as 


5(D) . , 

VFw2[C](D) 


( 2 ) 


where in the denominator the wavelet kernel is squared 
before performing the convolution. If the 7 -ray flux var¬ 
ied only on scales much larger than the extent of the 
wavelet kernel, and in the limit of a large number of 
photons, 5(D) would behave like a smoothed Gaussian 
random field. Consequentially, 5(D) can be loosely in¬ 
terpreted as the local significance for having a source at 
position D in units of standard deviations. 

As the wavelet kernel, we adopt the second member 
of the mexican hat wavelet family, which was shown to 
provide very good source discrimination power and 
which was used for the identification of compact sources 
in Planck data 0 - The wavelet can be obtained by 
a successive application of the Laplacian operator to a 
two-dimensional Gaussian distribution with width aiyR. 
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FIG. 1. SNR of the wavelet transform of 7 rays with energies 
in the range 1-4 GeV, 5(D). The black circles show the po¬ 
sition of wavelet peaks with 5 > 2; the red circles show the 
position of third Fermi-LAT catalog (3FGL) sources. In both 
cases, the circle area scales with the significance of the source 
detection in that energy range. The dashed lines indicate the 
regions that we use for the binned likelihood analysis, where 
latitudes |fe| < 2 ° are excluded because of the strong emis¬ 
sion from the Galactic disk. The subset of 3FGL sources that 
remains unmasked in our analysis is indicated by the green 
crosses. 


Here, Uh = 0.4° corresponds to the Fermi-LAT angu¬ 
lar resolution at 1-4 GeV, and F is a tuning parameter. 
We find best results when R varies linearly with latitude 
from R = 0.53 at 6 = 0° to F = 0.83 at 5 = ±12°. This 
compensates to some degree the increasing diffuse back¬ 
grounds towards the Galactic disk, while optimizing the 
source sensitivity at higher latitudes [doj - 

The resulting SNR of the wavelet transform 5(D) is 
shown in Fig. [TJ As expected, the Galactic diffuse emis¬ 
sion is almost completely filtered out by the wavelet 
transform, whereas bright sources lead to pronounced 
peaks. We adopt a simple algorithm for peak identifi¬ 
cation; we find all pixels in 5(D) with values larger than 
in the four adjacent pixels. We then clean these results 
from artifacts by forming clusters of peaks with cophe- 
netic distances less than 0.3°, and only keep the most 
significant peak in each cluster. 

In Fig. [ll we show the identified wavelet peaks with 
peak significance 5 > 2, as well as all 3FGL sources for 
comparison [l|. For sources that are bright enough in 
the adopted energy range, we find a good correspondence 
between wavelet peaks and the 3FGL, both in terms of 
position and significance (we compare the significance of 
wavelet peaks 5 with the 1-3 GeV detection significance 
for sources). 

It is worth emphasizing that for the adopted spheri¬ 
cally symmetric and centrally peaked distribution of the 
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CSP, most of the sources would be detected not directly 
at the GC but a few degrees away from the Galactic disk. 
This is simply due to the much weaker diffuse emission 
at higher latitudes. Our focus on latitudes |6| > 2°, thus, 
avoids regions where source detection becomes less effi¬ 
cient, due to strong diffuse foregrounds, without signifi¬ 
cant sensitivity loss for the source population of interest. 


3FGL sources. Before studying the statistics of the 
wavelet peaks in detail below, we remove almost all peaks 
that correspond to the known 3FGL sources based on a 
0.3° (1° for y/TS > 50) proximity cut. However, in order 
to mitigate a potential bias on Tmax, we do not mask 
peaks that correspond to 3FGL sources that are likely 
part of the GSP. We identify such MSP candidate sources 
by requiring that they (i) are tagged as unassociated, 
(ii) show no indication for variability and (iii) have a 
spectrum compatible with MSPs. The last criterion is 
tested by performing a fit of the above MSP reference 
spectrum to the spectrum given in the 3FGL (0.1-100 
GeV; five energy bins). Only the normalization is left 
free to vary. We require a fit quality of x^/DOF < 1.22 
(with DOF=4), corresponding to a p value > 0.3. 

We find 13 3FGL sources in the Inner Galaxy ROl 
that pass the above MSP cuts (listed in the Supplemen¬ 
tal Material). Interestingly, the average number of MSP 
candidate sources in same-sized control regions along the 
Galactic disk in the range \(.\ = 12°-60° is significantly 
smaller, with an average of 3.1. It is tempting to inter¬ 
pret this excess of MSP candidate sources in the Inner 
Galaxy as being caused by the brightest sources of the 
CSP, above the less-pronounced thick-disk population of 
MSPs 41,|^|. However, we emphasize that the status 
of these 13 sources is currently neither clear nor quali¬ 
tatively decisive for our results. Whether we mask them 
plays a minor role in the detection of the CSP below (but 
it does affect the inferred values for Tmaxl see Supplemen¬ 
tal Material). 


Statistical analysis. In Fig. [2] we show a histogram 
of the wavelet peaks in our ROl. We bin the peaks in 
a two-dimensional grid, which spans the projected angle 
from the Galactic Center 2°-17° and wavelet peak signif¬ 
icances in the range 1-10. The bin edges are as indicated 
in the figure. As expected, photon shot noise gives rise 
to a large number of peaks with low signihcances 5 < 3, 
and only a small number of peaks has 5 > 5. 

We assume that the number of peaks in each bin in 
Fig. [2] follows - in repeated experiments and random re¬ 
alizations of the CSP - to a good approximation a Pois¬ 
son distribution. We estimate the corresponding average 
number of expected wavelet peaks in each bin using a 
large number of Monte Carlo simulations, where we simu¬ 
late the diffuse background emission, random realizations 
of the source population and photon shot noise. 

In order to quantify what CSP luminosity function 
reproduces best the observations, we perform a binned 



Spatial and significance bins 


FIG. 2. Histogram of observed peaks in <S(n), in bins of 
projected radial distance from the GC and SNR values (black 
dots with statistical error bars). We show the expectation 
value for the case of a negligible CSP as blue bars, whereas 
the expectation values for the best fit are shown in red. 


Poisson likelihood analysis of the wavelet peak distribu¬ 
tion. The likelihood function is given by 

Ur ris 

^=nn l/iij (Z/niax, d^s)) ; (3) 

i=li=l 

where and rig are, respectively, the numbers of radial 
and peak SNR bins, Cij is the observed and /ty the ex¬ 
pected number of peaks, and P is a Poisson distribution. 
The expectation values depend directly on the maximal 
luminosity, Tmaxj as well as on the number of simulated 
sources, n. To ease comparison with the literature, we 
determine n as a function of $ 5 , which denotes the mean 
differential intensity of the CSP at 5 = ±5°, £ = 0°, and 
2 GeV. In the case of the GCE, this value was found to 
be i 1-5) X 10 “^ GeV“^ cm“^ s“^ sr“^ at 

95.4% C.L. [ 3 . 

Results. In Fig. we show the expectation values 
that we obtain when neglecting contributions from the 
CSP (and any other nondiffuse emission). This corre¬ 
sponds to good approximation to the case where the GCE 
is of truly diffuse origin, including the case of DM anni¬ 
hilation or outburst events. We find that the observed 
number of wavelet peaks with 5 < 2 is significantly lower 
than expected, whereas the observed number of peaks 
with iS > 3 is significantly higher. As we will show next, 
this is precisely the effect that is caused by a dim source 
population. 

We now turn to the case with a nonzero CSP con¬ 
tribution. In Fig. [31 we show the limits that we ob¬ 
tain on the two CSP parameters when fitting the his¬ 
togram in Fig. 13] as described above. We find that a 
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FIG. 3. Constraints on the maximum 7 -ray luminosity of the 
CSP, I/max, and the population averaged intensity at 6 = ±5°, 
I = Q°, and 2 GeV, 4 > 5 , as derived from our wavelet analysis. 
We show 68.7%, 95.4%, and 99.7% C.L. contours. We also 
indicate the values of "Fs where the source population can 
explain 100% of the GCE (horizontal gray band, 95.4% CL), 
and as vertical orange lines the luminosity of the brightest 
observed nearby MSPs. 


nonzero contribution from the CSP is favored at the 
level of at least 10 .Oct (when quoting the statistical sig¬ 
nificance, we conservatively take into account bins with 
5 < 5 only, which are most affected by a dim source 
population, and least affected by the masking of 3FGL 
sources [see the Supplemental Material for details]). The 
best-fit value for the total differential intensity is $5 = 
(9.0 ± 1.9) X 10“^ GeV“^cm~^s“^sr“^ and for the max¬ 
imum luminosity Lmax = (7.0 ± 1.0) x 10^^ ergs“^. As 
can be seen in Fig. [21 we obtain in this case a very good 
fit to the data. 

Our preferred range of the maximum 7 -ray luminosi¬ 
ties reaches up to Lmax < 1-04 x 10^^ ergs“^ (at 95.4% 
G.L), which is compatible with observations of nearby 
MSPs. We illustrate this by showing in Fig. [3] the 7 -ray 
luminosity of the brightest individually observed nearby 
MSPs as given in Ref. [s^ (we only show objects where 
second 7 -ray pulsar catalog [s^ distances are available; 
see Ref. [ 2 ^ for a detailed discussion about distance un¬ 
certainties). Furthermore, for the adopted slope of the 
luminosity function, a = 1.5, the best-fit value for the 
total differential intensity of the GSP $5 is consistent 
with the CSP accounting for 100% of the GCE emission. 

Discussion and conclusions. We found corroborating 
evidence for the hypothesis that the GCE is caused by 
a hitherto undetected population of MSP-like sources. 
We performed a wavelet transform of the 7 -ray emission 
from the Inner Galaxy, which removes Galactic diffuse 
emission and enhances point sources, and we studied the 


statistics of the peaks in this transform. We detected 
with 10 .Oct significance a suppression (enhancement) of 
low- (high-) significance wavelet peaks, relative to the 
expectations for purely diffuse emission. We showed that 
this effect is caused by the presence of a large number 
of dim point sources. The spatial distribution of wavelet 
peaks in the Inner Galaxy is compatible with a centrally 
peaked source distribution, and the inferred cutoff of the 
7 -ray luminosity function of these sources agrees with the 
observation of nearby MSPs. This source population can, 
for reasonable slopes of the luminosity function, account 
for 100% of the GCE emission. 

For the purpose of this Letter, which introduces a novel 
technique, we kept our analysis as simple as possible. 
In general, one might worry that our results could be 
affected by the presence of extragalactic and Galactic 
sources, by the thick-disk population of MSPs and young 
pulsars, by the details of masking and unmasking 3FGL 
sources, by the details of the adopted 7 -ray luminosity 
function, and by unmodeled substructure in the Galac¬ 
tic diffuse emission that is not removed by the wavelet 
transform. We address all of these points in the Supple¬ 
mental Material and show that it is rather unlikely that 
they affect our results qualitatively, although quantita¬ 
tive changes in the obtained best-fit values for $5 and 
Lmax are possible. In particular, we showed that the 
wavelet signal expected from the thick-disk population 
of MSPs is an order of magnitude weaker than what we 
actually observed and that interpretations related to un¬ 
modeled gas remain on closer inspection unlikely. 

The prospects for fully establishing the MSP interpre¬ 
tation within the coming decade are very good. Our re¬ 
sults suggested that upcoming 7 -ray observations with 
improved angular resolution (planned or proposed 7 - 
ray satellites like GAMMA-400 [i^, ASTROGAM, and 
PANGU 0 ) will allow us to detect many more of the 
bulge sources and study their distribution and spectra. 
For current radio instruments, it remains rather chal¬ 
lenging to detect a MSP population in the bulge (25| . 
but prospects for next-generation instruments are good. 
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Note added in proo/.-Recently, we became aware of 
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another group studying dim 7 -ray sources in the Inner 
Galaxy, using non-Poissonian photon statistics [i^ . 
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SUPPLEMENTAL MATERIAL 


In the Supplemental Material we discuss the possible impact of various systematic effects on our results. This 
includes a control region analysis, a discussion of various types of 7 -ray sources, substructure in diffuse emission, a 
thick-disk population of MSPs, sphericity and the role of negative wavelet peaks. 


A. Null results in control regions 


In order to estimate the effect of various systematic uncertainties, it is useful to apply our analysis on control 
regions along the Galactic disk (in the case of Galactic diffuse emission this was first systematically done in Ref. 0). 
Potentially unresolved substructure in the Galactic diffuse emission (e.g. in the form of giant molecular clouds, see 
below), and contributions from various Galactic and extragalactic source populations could be responsible for the 
detected wavelet signal in the inner Galaxy, but would in general also affect other regions in the Galactic disk. To 
this end, we focus on (partially overlapping) control regions along the Galactic disk, which are of the same size as the 
inner Galaxy ROI, but displaced by A£ = ±fc 20° with fc = 1, 2,3,4. 



Galactic longitude [deg] 


FIG. S-1. We show, as function of the central longitudinal position of the control regions (for ^ = 0° the main ROI), the 
significance of a CSP detection. We assume that the CSP is centered in each of the ROIs, and we refit the parameters Lmax 
and the number of sources. We indicate how different ranges of the wavelet peak significances contribute to the detection. All 
3FGL sources are masked in this plot, in order to be conservative. 


In Fig. IS-ll we show the TS value for a detection of the GSP for the main and the different control ROIs along the 
Galactic disk. We leave Lmax and t&s free to vary in each region independently. In the main ROI that covers the 
inner Galaxy we find the significant detection of a CSP that was discussed in the main text. As shown in the plot, 
this high significance is supported by the low, intermediate and high-significance SNR peaks of the wavelet transform 
separately. The directly adjacent regions also show relatively large TS values, which is either caused by the partial 
overlap of these control regions with the main ROI, or by a CSP that is more disk-like than assumed in our analysis. 
We will address the latter point below. However, in the outermost six control regions we find no significant detection 
of a CSP, for any of the considered values for Lmax and $5 (the large TS values at ^ = 80° are caused by one extremely 
bright source that generates fake peaks in its tails). This observation makes it already extremely unlikely that our 
findings are driven by a mismodelling of the local Galactic diffuse emission, or by extragalactic sources. We will 
address this in more detail below. 



























FIG. S-2. Similar to Fig. 3 in the main text, but showing limits derived for different ranges of the wavelet peak signihcances 
separately. The parameter degeneracies present in some of the fits cancel when fitting the entire significance range at once. We 
only show 68.7% and 95.4% CL constraints for clarity. 



FIG. S-3. Similar to Fig. 3 in the main text, but showing limits as derived for different spatial bins separately. 


B. Consistent wavelet signal in separate bins 

It is instructive to see how wavelet peaks with different significances contribute to the constraints on the luminosity 
function that we showed in Fig. [3] in the main text. To this end, we show in Fig. IS-21 the limits that we obtain 
separately from peak significances in the range S = 1-3, S = 3-5 and S = 5-10, respectively. All three constraints 
are mutually consistent to within Icr, leading to a consistent interpretation of the peaks shown in Fig. [T]in the main 
text. In all cases we find some degeneracy in the (Lmax, ‘I’s) plane. High significance peaks predominantly provide a 
stringent upper limit on Tmax, whereas low significance peaks mostly constrain the overall luminosity of the modelled 
source population. 

To show that our assumption on the spatial distribution of the CSP is consistent with the data, we show Fig. IS-31 
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FIG. S-4. Similar to Fig. 3 in the main text. We show the 68.7% and 95.4% CL contours for luminosity functions with spectral 
indices of 1.2, 1.5 and 1.7. 


the result obtained for the five different spatial bins independently. Ring 1-5 correspond tore [i°,i + 3°] with 
1 = 2,5, 8,11,14, respectively. We find constraints on Lmax and $5 that are mostly consistent to within Icr. 

Finally, we checked that the identified wavelet peaks are symmetrically distributed in the north, south, east and 
west parts of onr main ROI. Only at 5" > 3 we find a slight (statistically not very significant) asymmetry with more 
peaks in the south, which might be caused by the somewhat stronger Galactic foregronnds in the north, which makes 
point sonrce detection in the north more challenging. 


C. Mild dependence on the MSP luminosity function 


Theoretically, a is not well constrained and can plausibly range from 1.5 to 3, depending on the emission model [32 
Actnal MSP observations actually seem to indicate somewhat smaller values closer to a ~ 1.2 [s^. We show 
in Fig. lS-4l the 68.7% and 95.4% CL contours for different luminosity fnnctions, respectively with spectral indices of 1.2 
and 1.7. For a fixed intensity ( 4 * 5 ), hardening (softening) of the Inminosity function corresponds to an enhancement 
(suppression) of the number of sub-threshold point sources, which explains the direction in which the best fit region 
moves. We note that we obtain very similar TS values for all slopes that we considered. 


D. The role of unmasked 3FGL sources 

In the present analysis, we make use of the 3FGL, the third Fermi source catalogue, which is based on the first 
four years of Fermi pass 7 data. One important ingredient in onr analysis is the masking of 3FGL sources. These 
sources are of Galactic and extragalactic origin and leaving them unmasked would inevitably induce sizeable signals 
in our search for a sub-threshold source population in the bulge. However, as discussed in the main text, we keep 
nnassociated sources with MSP-like spectra unmasked. These sources could be part of the bulge MSP population 
that we are looking for, and masking them would bias our results. The 13 sources that pass our MSP cuts are listed 
in Tab. HI 

In general, falsely masking unassociated sources that actually belong to the bulge MSP population would push 
Amax to lower values, whereas falsely unmasking foreground sources would push it to large values. This is illustrated 
in Fig. IS-51 We show the case where we mask all unassociated sources, as well as the case where we adopt a weaker 
criterion for the spectral fit, leaving around 20 sources unmasked. We find that in both cases the best-fit value for 
Lmax moves in the expected direction, but the results remain consistent to within Icr. Furthermore, the significance 
of our wavelet detection that we qnote in the main text (where we inclnde 5 < 5 bins only) changes to 9.2cr when we 
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3FGL Name 


bn 

x7dof Trs 5 

L [lO®'^ erg/s] 

J1649.6-3007 

-7.99 

9.27 

1.07 

5.6 

7.4 

7.811'® 

J1703.6-2850 

-5.08 

7.65 

0.48 

2.4 

5.0 

3.47:1 

J1740.5-2642 

1.30 

2.12 

0.37 

6.4 

3.5 

14.917® 

J1740.8-1933 

7.43 

5.83 

0.77 

1.9 

2.3 

3-817 

J1744.8-1557 

11.03 

6.88 

0.40 

3.7 

3.4 

5.67® 

J1758.8-4108 

-9.21 

-8.48 

0.90 

5.6 

3.8 

4 . 87:1 

J1759.2-3848 

-7.11 

-7.43 

0.35 

4.6 

5.6 

5.97:1 

J1808.3-3357 

-1.94 

-6.71 

0.40 

6.9 

6.3 

8.07:1 

J1808.4-3519 

-3.15 

-7.36 

0.41 

4.6 

4.4 

5.07:® 

J1808.4-3703 

-4.68 

-8.19 

0.22 

4.9 

5.3 

4 - 37.1 

J1820.4-3217 

0.74 

-8.17 

1.04 

5.7 

1.7 

7.27:1 

J1830.8-3136 

2.35 

-9.84 

0.54 

5.9 

6.0 

5 . 07:1 

J1837.3-2403 

9.85 

-7.81 

0.28 

4.0 

3.0 

A ^+3.3 
‘ - 1.6 


TABLE I. List of the 13 unassociated 3FGL sources with MSP-like spectra, which we leave unmasked in our analysis. If the 
GeV excess is caused by dim point sources, it is likely that some or most of them are part of the GSP. The last four columns 
show the goodness-of-fit of the reference MSP spectrum, the 3FGL significance in the 1-3 GeV band, the corresponding peak 
of the wavelet SNR, and the y-ray luminosity (assuming 8.5 ± 2 kpc distance from the source and the reference stacked MSP 
spectrum from the main paper with a normalization that is obtained from a £t to the measured source flux). 



FIG. S-5. Similar to Fig. 3 in the main text. We show the effect of masking all 3FGL sources, or leaving a slightly larger 
number unmasked. 


mask all sources, and to lO.Scr when keeping 20 sources unmasked. The (un-)masking of 3FGL is hence not decisive 
for our qualitative findings, although quantitative results can be affected. 

It is interesting to note that the faintest 3FGL source in the inner Galaxy ROI that passes our MSP-spectrum cut, 
has a luminosity of L = 3.4 x lO^^ergs”^ if placed at 8.5 kpc distance. This is a good, though rough, indication for 
the de facto sensitivity threshold of Fermi-LAT for the detection of sources with MSP-like spectra in the bulge region. 

Lastly, one can use the sources in Tab.Uto compare the sensitivity of the 3FGL with our wavelet analysis. Averaging 
over the 13 sources, we find a ratio of S/\/TS ~ 1.0 ± 0.4, indicating that the sensitivity of the wavelet method is 
similar to the 3FGL sensitivity in a comparable energy range. Flowever, the scatter exceeds the one expected from 
statistical fluctuations alone, which can be attributed to differences in the systematics that affect the 3FGL and the 
wavelet analysis. 
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FIG. S-6. Latitude profile of sources, for different ROIs along the Galactic disk, as function of their longitudinal center (dotted 
lines). The ROIs have a size of 24° x 24°, with the Galactic disk, |fe| < 2°, masked. The different colors correspond to different 
source categories. In the case of unassociated sources, the solid (dashed) lines show only sources that pass the spectrum cut 
for MSP-like (young pulsar-like) sources. 
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FIG. S-7. Stacked histogram of wavelet peak values S that correspond to 3FGL sources in the inner Galaxy ROI. We show the 
contribution from different source categories separately. Unassociated sources generate mostly low-significance wavelet peaks, 
whereas for example source that are marked as pulsar in the 3FGL only generate peaks with a high significance (see discussion). 


E. The role of various source populations 

In Fig. IS-6[ we show the number of 3FGL sources in our inner Galaxy ROI as well as in various same-sized ROIs 
that are displaced along the Galactic disk in steps of A£ = ±24°. We show (identified and associated) extragalactic 
sources, various Galactic source classes and unassociated sources classes separately. We also show for the unassociated 
sourced how many sources pass our MSP cut. It is apparent that the number of unassociated sources strongly peaks 
in the inner Galaxy ROI, however with a clear asymmetry towards negative values of i, and another peak around 
i Ri 100°. However, after applying the MSP cut, predominantly sources in the inner Galaxy survive. In Fig. IS-71 on 
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the other hand, we show a histogram of the wavelet peak significance that our analysis attributes to the 3FGL sources 
in the inner Galaxy ROI. Again, unassociated sources play a major role and produce wavelet peaks down to values of 
5 ~ 1. On the other hand, 3FGL sources that are identified as pulsars appear only with 5 > 5 in our analysis. We 
will discuss in the following the potential impact of each of the source classes separately. 

Extragalactic sources. As shown in Fig. IS-61 the number of extragalactic sources in the individual ROIs along the 
Galactic disk fluctuates around values of about ~ 13. No significant suppression is observed in the inner Galaxy ROI, 
which would have indicated that it is more challenging to identify or associate extragalactic sources in this region. 
We will use here a simple argument to show that extragalactic sources cannot play a significant role for our results. 
The average number of 5 = 3-5 wavelet peaks in the above control regions along the Galactic disk is 20; in the inner 
Galaxy ROI it is 42 (we exclude here all 3FGL sources to be conservative). If extragalactic sources were the main 
contribution to these peaks, in addition to the about 10 wavelet peaks that are expected from statistical fluctuations 
alone (see Fig. [2] in the main text), the 42 observed peaks would constitute a > 5 (t upward fluctuation above the 
expected 20. This makes it extremely unlikely that extragalactic sources contribute significantly to our results in the 
inner Galaxy. Similar arguments can be made using wavelet peaks in the range S = 1-2. 

Supernova remnants and pulsar wind nebulae. As apparent in Fig. IS-61 only a very small number of sources along 
the Galactic disk at latitudes | 6 | > 2° (and almost none at | 6 | > 5°) are identified with supernova remnants or pulsar 
wind nebulae. In our ROIs their number is much less than the number of extragalactic sources, and their distribution 
is centrally not peaked, indicating that sources at these latitudes are mostly local. Sources in this category are 
typically more easily detectable at higher and lower energies than the energy range used in our analysis, and would 
be most likely listed in the 3FGL and hence masked if they were abundant and significant. We consider it hence as 
extremely unlikely that sources of this category significantly affect our results in the inner Galaxy. 

Young and millisecond pulsars. An interesting feature of pulsars in the 3FGL is that they always induce large 
wavelet signals in our analysis, as shown in Fig. lS-7| This makes indeed sense, since the identification of a 7 -ray source 
as pulsar requires the measurement of its pulsation, and hence a large enough number of photons. Furthermore, the 
pulsar energy spectrum often peaks close the energy range of our analysis. Already from Fig. IS-71 it is obvious that 
a large fraction of the unassociated sources, which mostly appear with lower significance peaks, must in fact be 
pulsars. Interestingly, as shown in Fig. IS-61 identified pulsars do not show a centrally peaked distribution, whereas 
the unassociated sources clearly do. This behaviour is expected, given that with increasing distance the pulsar 
identification becomes more challenging. 

Globular clusters. The 7 -ray emission from globular clusters that is not masked in our analysis, either because 
the globular clusters did not enter the 3FGL, or because they happen to be among our 13 unmasked unassociated 
sources, could in principle contribute to the detected signal. Since their total emission is usually due to several MSPs 
which appear as a single source for Fermi-LAT, their presence could bias Tmax towards larger values. However, given 
the simulation results from Ref. [ 2 ^ , we expect this effect to be small, and leave a more detailed discussion to future 
work. 

Unassociated sources. The peak of unassociated sources in the inner Galaxy, as shown in Fig. IS-61 appears clearly 
asymmetric, with a second peak at £ ~ 100°. As discussed above, a large fraction of these sources is expected to 
be young or millisecond pulsars. Obviously, the 3FGL unassociated sources do not directly contribute to our results 
since they are masked (except the 13 MSP candidates). However, the unassociated sources are extremely abundant 
even down to 5 ^ 1 , and their probable nature can be used as an indicator for what source population dominates 
just below threshold. 

In Fig. IS-61 we see that our MSP cut removes most of the unassociated sources, but leaves an excess of 13 unassoci¬ 
ated sources in the inner Galaxy, as discussed above. What is more, if we slightly modify the spectral criterion, using 
dN/dE oc e-^/ 2 GeV^- 2 ^ which is somewhat more pulsar-like (softer index, lower cutoff), this behaviour changes and 
we instead find excesses that are more correlated with the peaks of the unassociated sources away from the inner 
Galaxy. Although the statistical significance of this finding is rather difficult to quantify without a detailed study 
(which we leave for future work), this result is indicative. It suggests that a large fraction of the inner Galaxy unas¬ 
sociated (and sub-threshold) sources are likely MSPs, whereas unassociated sources in other parts of the disk have a 
larger fraction of young pulsars. The latter point is further supported by the fact that similar structures can be found 
in the longitudinal distribution of identified pulsars. 

In summary, we expect that our wavelet signal is dominated by whatever source class is responsible for most of the 
unassociated sources towards the inner Galaxy. Very likely, these are millisecond and young pulsars, with a somewhat 
higher MSP/young pulsar ratio than in the rest of the disk. Since these sources appear in general both in the Galactic 
disk as well as in the bulge, it is important to study whether the excess/suppression of wavelet peaks in the inner 
Galaxy points to a disk population, a bulge population, or to a combination of both. This will be discussed in the 
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£, Gal. longitude [deg] 


FIG. S-8. Similar to Fig. 1 in the main text (note the different color scale), but showing the transform of the diffuse BG model 
only, without Poisson noise being applied. Outside of the Galactic disk, |6| > 2°, which we exclude from our analysis, the 
significance of wavelet peaks remains below 0.5. The variance is below 0.1, which shows that even la peaks in the wavelet 
transform are unlikely to be strongly affected by the Galactic diffuse emission. 


section G below. 


F. Possible caveats concerning the Galactic diffuse emission 


In our Monte Carlo studies, we use the standard Fermi diffuse model for pass 8 data analysis. The wavelet 
transform of this model, without applying Poisson noise, but using the same exposure as in our main analysis, is 
shown in Fig. IS-81 Outside of the masked Galactic disk at | 6 | > 2°, we do not find any excesses with a significance 
larger than about 0.5cr. The main effect of such variations would be to offset the significance of random statistical 
fluctuations and sub-threshold point sources towards higher or lower values. This would not significantly affect peaks 
with a large SNR. However, it can potentially be important for low-significance peaks, because the collective shift of 
a large number of peaks, even by a small amount, could become statistically relevant. But since the variance of the 
SNR values shown in Fig. IS-81 is below 0.1, we do not expect that the details of the modeling of diffuse emission when 
doing MCs is going to affect our results. 

The Fermi diffuse models might not actually contain all relevant small-scale gas structures, and the effect of these 
missing structures on our results is not straightforward to estimate without a detailed analysis and modeling of the 
power-spectrum of gas at small scales. It is hence rather important that our non-detection of strong wavelet signals 
along the Galactic disk in Fig. IS-II largely excludes that mismodeling of local gas is the cause for the detected signal 
towards the inner Galaxy, since it would affect other parts of the disk as well. This is in particular true since there is 
relatively little molecular gas in our main ROI, compared to the control regions [i^. Thus, gas-related effects should 
be larger in the control regions than in the main ROI. 

If one insists on a gas-related interpretation, our results hence suggest that the wavelet signals are caused by 
unmodeled gas in the Galactic bulge, at a height of 0.3-1.5 kpc. If we assume a cosmic-ray density in the bulge 

s“^GeV“^ per hydrogen 


similar to the local one 
atom 
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the differential 7 -ray emissivity at 1 GeV is around 3 x 
47j . This implies that dense gas clouds with masses around 3 x 10^ Mq would be at 1 GeV roughly as bright 


as MSPs with a luminosity of T = 7 x 10^"^ ergs“^.[^ Interestingly, giant molecular clouds are known objects of 
that mass, and they can be dense enough to appear at GC distance point-like for Fermi. However, the scale height 
of known giant molecular clouds is at the level of a few 10 pc (they usually intersect with the Galactic disk) instead 
of the required ^ 1 kpc [i^. Furthermore, clouds of that size should give rise to CO emission in the range 0(10- 
100) Kkms“^, which is not seen in current observations ^j. If observed, such CO emission should be distributed 
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Galactic longitude [deg] 


FIG. S-9. Similar to Fig. IS-11 but including a thick-disk MSP population calibrated to local bright high-latitude MSPs as 
additional background. Only for illustration, we also show the effect of a 10x more dense thick-disk population (which is in 
contradiction with local observations). Note that we unmask 3FGL as described in the main analysis when deriving the wavelet 
peaks. 


north/south symmetric, as our wavelet peaks are too. 

If one could show that a large number of such giant molecular clouds (or other structures with similar mass and 
density) can form and be transported to kpc heights in the Galactic bulge, while hiding from all observations, the 
interpretation of the identified wavelet peaks in terms of unmodeled gas would remain a possibility. However, as 
of now, and for all of the above reasons, we regard gas-related interpretations of our results as rather unlikely and 
speculative. 


G. Potential impact of thick-disk population 


As argued above, the most relevant Galactic background in our ROI are expected to be pulsars, and in particular 
the MSP thick-disk population that reaches up to high latitudes. We will now show that a thick-disk population of 
MSPs (or other sources with a similar luminosity function) cannot be responsible for the observed signal. 

In most cases, the thick-disk population of MSPs is modeled as a cylindrically symmetric exponential distribution, 
with a scale height in the range 0.5-1 kpc and a scale radius of a few kpc, which is only poorly constrained by 
data (see e.g. Ref. [4l|). We will adopt here a distribution with a scale height of 1 kpc and a scale radius of 5 kpc, 
which was previously used to argue against the MSP-origin of the Fermi GeV excess [sH. The distribution reads 
n oc exp {—R/Rs) exp (— jzj /Zg), with i?s = 5 kpc and Zs = 1 kpc. We will address below how the results change when 
other parametrizations are adopted. 

As 7 -ray luminosity function, we adopt an inverse power-law with Lmin = 10^^ ergs“^, Lmax = 7 x 10^^ ergs“^ and 
index a = 1.5. We fix the overall normalization of the disk source density such that the number of bright MSPs at 
high latitudes, |6| > 15°, is consistent with the number of such MSPs listed in the 3FGL. As flux threshold for bright 
MSPs we adopt a flux that corresponds to a y-ray luminosity of 10^^ ergs“^ at 3 kpc distance (9.2 x 10“^^ ergs“^cm“^ 
in our energy range). We find 31 MSPs above that threshold flux and note that since the number of unassociated 
high-latitude bright nonvariable sources with a curved enough spectrum is small, this number cannot increase by more 
than 50% when more unassociated sources are identified as MSPs) [l[. For the present scenario, we find that the total 
number of thick-disk sources with y-ray luminosity above 10^^ ergs“^ is ~ 30000. 

Within 2 kpc of the Galactic center, this thick-disk population predicts around 1300 MSPs, which is more than an 
order of magnitude below the number that we find in the best-fit scenario for the bulge population (around 35000 
MSPs above 10^^ ergs“^). This implies, as already argued in Ref. 3l|, that a thick-disk population with the adopted 
geometry cannot be responsible for the Fermi GeV excess. However, it also trivially implies that the number of 
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wavelet peaks caused by thick-disk sources in the inner 2 kpc is about an order of magnitude below what is predicted 
by our best-fit bulge population, and hence an order of magnitude below what is actually observed. This still leaves 
the possibility that thick-disk MSPs on the line-of-sight towards the inner Galaxy, outside of the inner ^ 2 kpc, could 
affect our results. We will discuss this next. 

Within our ROI, the thick-disk population predicts 3.3 sources outside of the inner 2 kpc with a flux in the range 
(4.6-7.7) X 10“^^ ergs“^cm“^ (this corresponds roughly to (4-7) x lO^^ergs”^ when the sources are put at 8.5 kpc 
distance). These sources could reasonably contribute to wavelet peaks in the 3 < 5 < 5 range. The actually observed 
number of peaks in that range above the null hypothesis is about 35. It is hence clear that foreground sources from 
the above thick-disk population cannot cause the observed signal. One might think of two ways around. 

First, one could reduce the scale radius of the thick-disk population such that the number of sources in the inner 
2 kpc increases by a factor around ten (scale radii around 1-2 kpc could do the job). This would give rise to a 
wavelet signal similar to what is observed. However, such a population would also predict a significant diffuse 7 - 
ray emission similar to the level of the Fermi GeV excess, just with a morphology that is incompatible with the 
observations. A population with a scale radius of 1-2 kpc would indeed commonly be referred to as bulge population. 
Such a population would be very similar to the bulge population that we put forward in the main part of the paper 
as explanation for the Fermi GeV excess, with the main difference being that our population fits better the excess 
morphology. 

Second, one could increase the number of MSPs in a ring-like region around the Galactic bulge, excluding the inner 
2 kpc, such that these additional ring-like distributed sources will enhance the number of foreground sources without 
affecting the number of sources in the Galactic bulge. In this case, however, the wavelet signal should clearly be more 
extended along the Galactic disk than what is shown in Fig. IS-11 since such a ring would not be centrally peaked 
and extend to longitudes of at least ~ 25°. For illustration, we here quote the relative number of wavelet peaks one 
expects in the control regions along the disk and the main ROI produces by such a ring (1 kpc scale height, 5 kpc 
scale radius, the inner 2 kpc radius excluded): A£ = {±80, ±60, ±40, ±20,0} and iVpeaks oc (1.7, 2.4,3.5,4.9,3.6}. 
Moreover, in order to avoid a conflict with the above calibration with bright high-latitude sources, the ring should 
be further constrained to lie within < 5 kpc, which however would still leave a too flat central distribution of wavelet 
peaks. 

Finally, we show in Fiar. lS-Ql how the TS values are affected if the thick-disk population (or a lOx denser population) 
is added as an additional background component. For simplicity, we assume that the thick-disk population causes 
deviations of the expectation values fiij in Eq. ([3|) from the null hypothesis that are proportional to the deviations 
caused by the best-fit bulge population. We adjust the normalization of these deviations such that the number of 
additionally predicted 3 < 5 < 5 peaks in the main ROI is 3.3 (as motivated by the above discussion). We then 
add these thick-disk-induced deviations from the null hypothesis as negative and positive contributions to the model 
predictions in Eq. and repeat the CSP fit to inner Galaxy data. We repeat this procedure in all of the control 
ROIs used in Eig. IS-ll reweighing the thick-disk contribution properly at different Galactic longitudes. Erom Eig. IS-91 
it is clear that only a thick disk ten times denser than what is actually observed at higher latitudes could significantly 
affect, although not completely remove, the excess of wavelet peaks in the inner Galaxy. 


H. Further discussions 

In order to test whether the spatial distribution of wavelet peaks in our analysis is indeed compatible with a 
spherically symmetric distribution, we re-binned the wavelet peaks into a north/south region, defined by 12 ° > | 6 | > 
max(|£|,2°), and an east/west region, defined by 12° > \£\ > | 6 | > 2°. For each of the two regions, we derive the best- 
fit values for Lmax and $ 5 . The results are shown in Fisf. IS-lOl We find that the inferred parameters are consistent 
within one sigma, with a slightly stronger signal in the north/south region. This indicates that the excess of wavelet 
peaks, if interpreted in terms of a bulge source population, is consistent with a spherical distribution of these sources. 

Given that unresolved sources only add positively to the Galactic diffuse emission, whereas a mismodeling of the 
gas could cause both positive and negative variations, it is tempting to think that a detection of negative wavelet 
peaks would disfavour an interpretation in terms of unresolved point sources. In fact, we do find a suppression of 
—2 < iS < — I peaks, and an enhancement of —4 < S < —3 peaks in the inner ROI when searching for negative instead 
of positive peaks. But unfortunately, this cannot easily be used to discriminate diffuse modeling artefacts (which, as 
we discussed above, are anyway unlikely, as they should show up in the entire disk) from sub-threshold point sources. 

Maybe somewhat un-intuitively, negative wavelet peaks can indeed be generated by a large number of weak, positive 
point sources. This happens in the tails of our simulated sources, where the wavelet transform becomes negative (this 




16 



Maximum 7 -ray luminosity, I/max [erg s 


FIG. S-10. Similar to Fig. [3] in the main text, but derived in redefined regions within our main ROI that are useful to check 
for the sphericity of the signal (we still mask |&| < 2°). We find best-fit parameters for Lmax and $5 in the north/south and 
east/west regions that are consistent within one sigma. 


effect is visible as rings around bright sources in Fig. 1). We estimated the expected number of negative wavelet 
peaks for the best-fit scenario in Fig. 3 by Monte Carlo simulations, and find results that are completely consistent 
with the observed number of negative peaks. However, given that the number of positive and negative wavelet peaks 
are correlated (they are caused by the same sources), one cannot easily use observations of negative wavelet peaks to 
further constrain the model parameters. This, and the fact that an appropriate masking of 3FGL sources (including 
also the ring around each source) reduces significantly the effective size of the ROI, make an efficient use of negative 
peaks in our analysis difficult. However, it is re-ensuring that both the observed negative and positive wavelet peaks 
are consistent with the respective predicted number of negative and positive wavelet peaks for the same sub-threshold 
point source population. 

Finally, we briefly comment on the recent analysis of the Galactic center data by the Fermi-LAT collaboration 
[l 6 | and compare their use of wavelets to ours. Ref. [l^ uses wavelets to find seeds for the identification of point 
sources (also see Ref. [s^ for details about the adopted method), followed by standard maximum-likelihood fits for 
further source identification. This can be potentially affected by interstellar emission modelling. We instead study 
the statistics of local-maxima in the wavelet-transformed sky map, which is largely background-model independent 
(although small scale fluctuations could in principle be relevant, as discussed above). We find good correspondence 
between our wavelet peaks and the 3FGL catalogue (see Fig.[T|), which supports the validity of our approach. However, 
we also find that this agreement and the quality of the wavelet analysis in general critically depends on the ado pted 
wavelet type and size. All these points make it difficult to directly compare our results with those of Ref. [I 6 | . 
However, we note that the Galactic disk, where Ref. 0 finds that their identified sources most strongly trace the 
edges of the interstellar emission and thus might constitute false positives due to gas fluctuations (see their Fig. 8 ), 
is masked in our analysis. 






